@ ------------------------------------------------------ @
@ GAUSS code for aggregate income model with IID errors  @
@ ------------------------------------------------------ @

/* First: run stage 1 panel IID/ARCH/GARCH to put lambdat and gmm in memory */

@ ++++++++++++++++++++++++++++++++ Estimation +++++++++++++++++++++++++++++++++++++ @

par_agg=2;     
par_in1_agg=0|-1;
{xout,fout,gout,cout}=optmum(&li_agg,par_in1_agg);
par_final_agg=trans_agg(xout);
kap=par_final_agg[1,1];  
varpi=par_final_agg[2,1];

if cout==0;
LL_val_agg=-fout;
hessian_agg=hessp(&li_agg,xout);
if minc(eigrs(hessian_agg))<=0.00000001;
cls;
print "   ";
print "Numerical Hessian using 'hessp' is not positive definite, using '_opfhess' instead";
omega_hat_agg=inv(_opfhess);
else;
cls;
omega_hat_agg=inv(hessian_agg);
endif;
endif;

grdn_agg=gradp(&trans_agg,xout); 
cov_agg=grdn_agg*omega_hat_agg*grdn_agg';     
sd_agg=((sqrt(diag(cov_agg))));
AIC_agg=(-2*LL_val_agg+2*par_agg)/25;

@ +++++++++++++++++++++++++++++++ Output ++++++++++++++++++++++++++++++++++++++ @

gmm=fin_agg(par_final_agg);

sqngamma=sqrt(n)~sqrt(n);
par_gamma=kap|varpi;

@ ++++++++++++++++++++++++++++++++ Report +++++++++++++++++++++++++++++++++++++ @

print "  ";
print "********************************** Estimation Output *****************************************";
print " ";
print "Kappa and unconditional variance aggregate shock (left) and standard errors (right)";
print par_final_agg~sd_agg;
print "**********************************************************************************************";
print "  ";
print "Maximized log-likelihood (left) and Akaike information criterion (right)";
print LL_val_agg~AIC_agg;
print "  ";

@ ++++++++++++++++++++++++++++++++ Likelihood +++++++++++++++++++++++++++++++++++++ @

proc li_agg(prmtr);
local par_in,kap,varpi,sum,sumli;
local i,Y,vt,vt2;
 
par_in=trans_agg(prmtr);
kap=par_in[1,1];
varpi=par_in[2,1];

sum=0;
Y=lambdat;
vt=0;
vt2=varpi;

i=1;do until i>25;

vt=Y[i,1]-kap;
vt2=vt^2;

sum=sum+ln(varpi)+vt2*inv(varpi);

i=i+1;endo;

sumli=-(-0.5*(25*ln(2*pi))-0.5*sum);

retp(sumli);
endp;

@ ++++++++++++++++++++++++++++++++ Data +++++++++++++++++++++++++++++++++++++ @

proc fin_agg(prmtr);
local par_in,kap,varpi;
local i,Y,vt,vt2,shock_agg,variance_agg;
 
par_in=prmtr;
kap=par_in[1,1];
varpi=par_in[2,1];

Y=lambdat;
vt=0;
vt2=varpi;
shock_agg=zeros(25,1);
variance_agg=zeros(25,1);

i=1;do until i>25;

variance_agg[i,1]=varpi;
vt=Y[i,1]-kap;
shock_agg[i,1]=vt;
vt2=vt^2;

i=i+1;endo;

i=1;do until i>kkk;

year=gmm[i,2];
gmm[i,7]=shock_agg[year-67,1];
gmm[i,8]=variance_agg[year-67,1];

i=i+1;endo;

retp(gmm);
endp;

@ ++++++++++++++++++++++++++++++++ Transformation +++++++++++++++++++++++++++++++++++++ @

proc trans_agg(x0); 
local x1;

x1=x0;
x1[2,1]=exp(x0[2,1]);

retp(x1);
endp;

